Approximating the ground state of fermion system by multiple determinant states: 

matching pursuit approach 
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We present a simple and stable numerical method to approximate the ground state of a quantum 
many-body system by multiple determinant states. This method searches these determinant states 
one by one according to the matching pursuit algorithm. The first determinant state is identical 
to that of the Hartree-Fock theory. Calculations for two-dimensional Hubbard model serve as a 
demonstration. 

PACS numbers: 02.70.-c, 31.15.Ar, 71.15.-m 



Searching a single determinant state to approximate a 
quantum ground state, namely, the Hartree-Fock (HF) 
algorithm, plays important role in the understanding of 
nuclear, atomic, and molecular structures. It is a long 
standing effort to extend the HF theory into a truly 
first principle method by searching multiple determinant 
states to span a quantum state, for recent examples, 
see [H, 0, y, 0, HI and references therein. The attract- 
ing feature is that this approach is very stable and free 
from the sign problem. It in principle can apply to a wide 
variety of systems. However, first principle calculation in 
terms of multiple determinant states is still a challenge. 
In fact, including multiple determinant states in the vari- 
ational treatment, which is the common approach, often 
results in very complicated formulations. The compu- 
tation cost is usually impractically demanding. Some 
realistic implementations impose restrictions on the de- 
terminant states. For example, the Multi-configuration 
Hartree-Fock theory [l|, a time dependent extension to 
the HF theory, requires single particle states to be or- 
thogonal with each others. Here, we use a new approach 
to approximate the quantum ground state via multiple 
determinant states. 

Here, based on the matching pursuit (MP) algo- 
rithm [a, 0, H, Q , we show a numerical method to search 
determinant states to span the ground state of a fermion 
system. The determinant states are found one by one 
from all possible determinant states. Searching the first 
determinant state is identical to the Hartree-Fock theory. 
A significant feature of the current method is that several 
tens of basis determinant states are enough for reasonable 
result, and one can reach high accuracy by searching one 
or two thousands basis states. These numbers of basis 
states are several orders smaller than that of the Stochas- 
tic diagonalization algorithm which searches orthog- 
onal determinant basis states stochastically to span a 
quantum wave function. In comparison with other al- 
gorithms of search determinant states to span a fermion 
ground state, such as the Path-Integral Renormalization 
Group (PIRG) algorithm 0, [|, the current MP based 



method is quite simple and efficient. 

The MP algorithm is originally designed for signal pro- 
cessing [6] . It is now popular on the engineering commu- 
nity for coding, analysis, and compression of video and 
audio data @,H, @ • This algorithm searches some basis 
states from an over-complete basis set to represent a se- 
quence of data. The basis states are found one by one. 
The convergence of the MP algorithm is proved mathe- 
matically. For sufficient redundancy of the over-complete 
basis set, the convergence can be exponential [§]. The 
MP algorithm is insensible to the dimension of the data, 
and thus promises applications in quantum many-body 
systems. In Ref. the authors employ this algorithm 
to propagate quantum wave functions via split operator 
method in the Gaussian wave packet basis. A encourag- 
ing result is that several tens of Gaussian wave packets 
are able to accurately represent quantum wave function 
of a 20-dimensional model. 

The goal of MP algorithm is to obtain a sparse rep- 
resentation of a signal. To represent a quantum many- 
body wave function, ip, the MP algorithm searches an 
over-complete basis set and finds some basis states, <fii, 
(f>2, •■•,</>„, such that the combination of the basis states, 
ijj n = oi\<\>\ + ■ • • + a n <p n can best approach the state tp. 
Mathematically, that is to require ip n has minimum dis- 
tance with -0, i.e. \ip — ij) n \ reaches minimum. The basis 
states are found one by one. At /c-th step, the basis state 
4>k is obtained such that the combination of the basis 
states ipk = ol\<\>\ + • ■ • + ak4>k has minimum distance 
with the state ifi, i.e. \ip — ipk\ has minimum for all pos- 
sible choice of cj>k- Each more step brings the ipk closer 
to the target state ip, i.e., the distance \ip — ipk\ decreases 
with k. 

The eigenvalue problem is equivalent to find minimum 
values of the Rayleigh quotient 



E= MiJ|V)/W#>, 



(1) 
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where H is the Hamiltonian and %j) is the trial wave func- 
tion. Calculation of the ground state by the MP algo- 
rithm is to search some basis states to span the ground 
state. The basis states are found one by one from an over- 
complete basis set. Each searching process obtains one 
basis state such that the combination of this basis state 
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and those already found ones minimizes the Rayleigh 
quotient for all possible choice of the current basis state. 
This process of finding a new basis state continues until 
convergence of E. Without loss of generality, in the fol- 
lowing discussions, we focus on fermion systems, and use 
all possible Slater determinant states as over-complete 
basis set. 

Note that, for fermion systems, the first step is to find 
a Slater determinant state that minimizes the Rayleigh 
quotient. This is just the well known Hartree-Fock ap- 
proximation. We employ an iterative method to search 
a new determinant state, including the first one. We de- 
note the single particle basis states as (i = 1, • • • , n), 
and af (aj) the operator for creation (annihilation) of 
the state \i), i.e., \i) = af\0) with |0) the vacuum state. 
A determinant state can be expressed as 

m 

where m is particle number and Fj~ (Fj) is creation 
(annihilation) operator for single particle state, Fj = 

c ij a i + ' ' ' + c njO-n- Searching for the determinant state 
0) is equivalent to find the coefficients {cy} (or the op- 
erators {Fj'}). 

We use an iterative relaxation procedure to search the 
operators {Fj }■ From an initial trial state in the form of 
(|2"|) which can be chosen randomly, we optimize F^~ , F^~ , 
• • • , F+ consecutively. Each step of the optimization 
lowers the Rayleigh quotient. This iteration continues 
until the convergence of the Rayleigh quotient. Note that 
the determinant state © is a multi-linear function of 
the coefficients {c^}. For a fixed j, \<f>) is just a linear 
function of C\j , ■ ■ • , c n j : 

\4>) = ^djl&j), (3) 

i 

where \4>ij) — d\<p)/dcij. Thus an approximate ground 
state = J2i OLi(jft> + acj) can be written as 

= 2 + H ac iJ<l>ij ■ ( 4 ) 

i i 

This means that we can improve and hence up- 
date the operator Fj by finding the lowest eigenstate 
of the Hamiltonian in the subspacc spanned by {|<^y), 
i = 1, ••• ,n} and those previously found determinant 
states |</»W>. 

Such relaxation procedure to update the operators Fj 
is the key ingredient of this contribution. Suppose we 
have already obtained fc — 1 determinant states \4>^'}, the 
searching process for fc-th determinant state \4>^) = \4>) 
in the form involves the following iteration: 

(1) randomly generate a determinant state |</>). 

(2) For j = 1, 2, • • ■ , m, do the following iteration loop 
to update \(j>): 



(2a) Calculate the matrix elements of the Hamiltonian 
in the subspace sj fc ^ spanned by {l^' 1 '), • • • , \^ k ~^), 

\<t>ij)> ■■■ i l<?W)}; 

(2b) Find the ground state Vf^ of the Hamiltonian in 
the above subspace Hp, ^ = £\ + J2i Pij&j'i 

(2c) Update Fj by setting Cjj = 0^ (i = 1, • • • ,n); 
Then make ^ + |0) orthogonal to other single particle 
states {i^ + |0), I 5^ j}, and restore Fj |0) to unit length 
by a normalization procedure. 

(3) Check the convergence of the Rayleigh quotient. 
Repeat the step (2) until reaching convergence. 

In case of k = 1, the above searching process of finding 
the first determinant state is the same as the Hartree- 
Fock algorithm. A randomly generated initial trial state 
needs several tens of iteration rounds to converge. Each 
of subsequent determinant states needs about similar 
rounds of iteration. Here, the main numeric cost is the 
step (2a) for calculation of the matrix elements of the 
Hamiltonian between basis states. The step (2b) of find- 
ing lowest eigenstate in the subspace can be implemented 
efficiently via iteration algorithm [ll|, [12] that needs only 
small portion of the computation cost. 

Starting from k = 2, the number of iteration to ob- 
tain cj)k depends on the initial choice. If a trial state 
has large overlap with the state (H — Ek-i)Y& k-i) , one 
may reach convergence by just a few rounds of iteration. 

Here, Ek-i and ^k-i = Y^i=i a i4>^ are the approxi- 
mate ground state energy and wave function obtained 
in the previous step. One can understand this property 
by considering minimization of the Rayleigh quotient in 
the two dimensional subspace spanned by ^fe-i and the 
trial state <f> Q. From this observation, we perform a 
preparing treatment of the trial state before step (2) of 
the above iteration procedure. 

The preparing treatment of the initial trial state <p is 
to modify the state <fi so that it has maximum overlap 
with the state (H — Ek-i)\&k-i)- This procedure is easy 
to carry out by exploiting the fact that state (f>, or the 
overlap (4>\H — Ek-i^k-i)) is a multi-linear function of 
the coefficients (or the operators Fj~). We maximize 
the overlap iteratively by updating the operators Fj, 
(j = 1, • ■ • ,m), consecutively. Usually, 3 to 5 rounds 
of the iteration are enough. After such preparing treat- 
ment, one usually needs about 2 to 3 iteration rounds 
of the searching process to minimize the Rayleigh quo- 
tient. Thus, such preparing treatment makes the overall 
procedure about 5 to 10 times faster. 

As an optional choice to achieve high accuracy, one 
can perform backward optimization after reaching con- 
vergence in the above procedure. This procedure updates 
the already found basis states one by one (One can also 
choose to update some selected basis states Q). The op- 
eration to update a basis state is the same as searching 
a new basis state. It is numerically expansive to per- 
form the backward optimization. In fact, searching the 
basis states one by one is a kind of restriction on the de- 
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terminant states, and the backward optimization means 
removing such constraint. 

At first sight, the current method shares some features 
with the Path-Integral Renormalization Group (PIRG) 
algorithm 0, However, based on different strate- 
gies, the PIRG and the current method are two differ- 
ent methods of searching basis determinant states. The 
PIRG filters out the ground state by repeatedly expand- 
ing e~ TH \ip) into summation of determinant states and 
keeping some of the determinant states as new basis 
states to update the trial ground state \ip). At each step, 
the PIRG must update whole basis states, while the cur- 
rent method only adds (or updates) one basis state via re- 
laxation method and exploiting the multi- linearity of the 
determinant states. Updating basis states in PIRG, i.e., 
choosing some determinant states from those ones that 
span the state e~ rH \tp), involves diagonalization of many 
sizable matrices. The diagonalization in subspace is a 
major numeric cost of PIRG, while it takes only a small 
portion of numeric operations in current method. The 
current method only calculates matrix elements of the 
Hamiltonian between determinant states, this is much 
easier and more efficient than expanding e~ rH \i/j) (or 
H\ip)) into summation of determinant states. 

We test the above method via the two dimensional 
fermionic Hubbard mode on a N = L x L square lattice 
with periodic condition. The Hamiltonian reads 



(ij) 



3& 



i") + U^n^n^. (5) 



Here Cj a (cj a ) is the creation (annihilation) operator of 



an electron with spin a at j-th site and rij a = Cj a Cj a . U 
is the on-site Coulomb energy. The summation (ij) runs 
over nearest-neighbor sites. 



TABLE I: Ground state energies of some the 4x4 systems. 



system U/t CPMC SD 



MP 



N 



Exact 



10/16 
14/16 
16/16 
10/16 
14/16 
16/16 
14/16 
16/16 



4 

4 

4 

8 

8 

8 

12 

12 



-19.5808 
-15.7296 

-17.4800 
-11.648 

-9.696 



19.58 
15.49 
13.59 
17.40 



19.5775 
15.7107 
■13.5963 
■17.4625 
•11.6763 
•8.41565 
9.79500 
5.95238 



1874 
1865 
1622 
2000 
2000 
1850 
2000 
1950 



-19.5808 
■15.7446 
•13.6219 
-17.5104 
-11.8688 
-8.46889 
-10.0515 
■5.99222 



Table U shows ground state energies (in the unit of t) of 
some 4x4 systems. The column "system" indicates the 
number of electrons versus the lattice number. N is num- 
ber of basis determinant states to span the ground state 
of the current method (MP). We list the results of the 
Constrained Path Quantum Monte Carlo (CPMC) [ill ]. 
Stochastic Diagonalization (SD) Q, and exact diagonal- 
ization [IJ, [U [r| [TtJ for comparison. The accuracy 
of our method is almost unchanged for various interac- 
tion strength U and filling number. This demonstrates 



the stability of the method. All initial trial states for 
searching basis determinant states are randomly gener- 
ated without any symmetry consideration. We use con- 
vergence rate e = 2\E n — E n -i\/\E n + E n —\\ to deter- 
mine the number of basis states, where E n and E n —\ are 
ground state energies obtained with n and n — X basis 
states, respectively. The searching for basis states stops 
if e is smaller than a criteria eoj or maximum acceptable 
number of basis states is reached. Usually, eo = 1CP 5 is 
enough to obtain quite reasonable result. At this setting, 
one usually needs several hundreds basis states which 
increases slowly with U. If eo = 10~ 6 , one needs sev- 
eral thousands of basis states for convergence. Roughly 
speaking, result from about 100 basis states is quite well. 
The interest point is that the beginning several tens of de- 
terminant states, usually less than 60 basis states, make 
dominant contribution. And the first one, i.e., the HF 
approximation, contributes most. The number of dom- 
inant basis states increases slowly with the interaction 
strength U. After the dominant basis states, the contri- 
bution from each of following ones drops rapidly. 



TABLE II: Correlation functions of some the 4x4 systems. 



Method system U/t p(2, 1) 


S m (ir,n) 




CPMC 10/16 


4 


-0.0556 


0.731 


0.504 


MP(800) 10/16 


4 


-0.0559 


0.7315 


0.5089 


Exact 10/16 


4 


-0.0556 


0.73 


0.506 


QL 16/16 


4 


-0.0363 


3.42 


0.3946 


MP(437) 16/16 


4 


-0.0478 


3.490 


0.3888 


Exact 16/16 


4 


-0.0475 


3.64 


0.385 


CPMC 10/16 


8 


-0.0462 


0.761 


0.4403 


MP(800) 10/16 


8 


-0.0493 


0.7645 


0.4412 


Exact 10/16 


8 


-0.0485 


0.75 


0.443 



Table [H] shows ground state's correlation functions of 
some 4x4 systems. The comparing results of CPMC, 
Quantum Langevin (QL), and exact diagonalization are 
from [l3j], [l7[, and 3117 1, respectively. Here Sm and 
Sd are magnetic and density structure factors [13] , re- 
spectively; and p(r) is the one body density matrix. The 
number in the column "method" is the number of basis 
states of our method (MP). The current method obtains 
the ground state energy and wave function at the same 
time. Then calculation of the correlation functions and 
other related quantities is a trivial task that simply reads 
the wave function. In comparison with the exact result, 
we see that the wave functions are almost in the same ac- 
curacy as the correspondent energies. Again, the begin- 
ning several tens of basis states make major contribution. 
To demonstrate this property, as indicated in the paren- 
theses, we use only several hundreds of basis states to 
calculate the correlation functions in table [TXJ Since our 
program does not perform any symmetry treatment, we 
can only compare non-degenerated ground states with 
the exact result. In fact, the present method can take 
into account of symmetries. After the searching process 
of finding basis determinant states, the resultant wave 
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function is usually a combination of ground states with 
different symmetries. One may employ, e.g., the project 
technique [HI to filter out the target symmetry. This 
may further improve the accuracy. 



TABLE III: Ground state energies (in the unit of t) of some 
large systems with U = At. 



system 



QMC VMC SD MP 



N 



26/6 x 6 
34/6 x 6 
36/6 x 6 
50/8 x 8 
54/8 x 8 
62/8 x 8 
64/8 x 8 
100/10 x 10 



-42.32 
-33.30 
-30.96 
-72.80 
-67.55 
-57.70 
55.23 
-86.70 



■32.76 
■30.384 



-40.77 



-67.00 



41.0757 400 
32.7323 940 
-30.5166 905 
68.5029 630 
63.8981 560 
55.5255 490 
■53.583 472 
82.9549 152 



0.16 



0.12 



LU 



0.08 



0.04 



10 u 



10' 





16/16 U=12 — •— 




1 6/1 6 U= 4 -o- 




36/36 U= 4 -■- 




54/64 U= 4 - 




O A If? A 1 1 A a 

b4/b4 U- 4 — * — 




\ 100/100 U= 4 — •— 





n 



10" 



10 J 



FIG. 1: Relative error E r versus number of basis states n for 
various filling numbers and system sizes. 



Table Hill shows ground state energies for large system 
size that exact diagonalization is impossible. Here N is 
the number of basis states of the current method (MP). 
Our result is quite close to that of the SD [f|, and the 
variational quantum Monte-Carlo (VMC) 0]. It is worth 
to note that the number of basis states of our method is 
several orders smaller than that of the SD algorithm. As 
a consequence, our method needs much less memory, and 
there is no need for external storage. There are several 
percent of discrepancy with the Quantum Monte Carlo 
(QMC) result this disagreement increases with the 
system size. This needs further investigations. The dis- 
crepancy between QMC result and the strict variational 
result is also found by other authors, see, e.g., (3. H. [l8| . 
For practical applications, the extrapolation method in- 
troduced in PIRG's implementation is an useful tool to 
handle the discrepancy with QMC results [H, [l|| . Simi- 
lar to the 4x4 cases, the beginning several tens of basis 
states make dominant contribution. For a fixed interac- 
tion strength U, the contribution from the dominant ba- 
sis states increases with the system size. This means that 
one needs less basis states for larger system size. On the 
other hand, the computation cost to search a basis state 
scales about quadratically with the system size. The role 
of preparing step for searching a basis state is more signif- 
icant for larger system size. With out this preparing step, 
after several tens of dominant basis states, the overlap be- 
tween a randomly generated trial state <j> and the state 
(H — almost vanishes. One must substan- 

tially increase accuracy requirement for following steps, 
which in turn increases numeric cost. The computation 
cost scales about quadratically with the number of basis 
states. 

Fig. Q] shows the relative error E r = \(E n — E )/E \ 
versus the number of basis states n, where E n is the 
ground state energy obtained with n basis states, and E 
is the converged ground state energy (or exact ground 
state energy if available). This illustrates the overall 
properties of the method. The first basis state makes 
most important contribution. It accounts for about 80% 



to 95% of the ground state energy depending on sys- 
tem size and correlation strength. Roughly speaking, 
the mean field effect increases with the system size, and 
decreases with the correlation between electrons of the 
system. For system size A^o = 4 x 4 with U — 12t, the 
contribution of the first basis state is about 80%. As the 
system size reaching Nq = 10 x 10 with U = At, the con- 
tribution of the first basis state is more than 90%. As a 
consequence, for a same accuracy requirement, the num- 
ber of necessary basis states decreases with system size. 
The convergence rate is fast for the beginning several tens 
of dominant basis states. These basis states contribute 
more than 95% to the ground state for moderate correla- 
tion. Then contribution from each of the following basis 
states drops rapidly. However, the major computation 
cost for Fig. Q]is to search the remaining basis states. In 
practical calculations, one usually needs several tens of 
basis states for a reasonable accuracy. 

We perform backward optimization for some cases to 
improve the accuracy. There is very limited improve- 
ment from the backward optimization if the MP method 
is converged. The improvement is usually less than 0.5%. 
If backward optimization is performed before the conver- 
gence of MP process, the improvement can be more than 
1.5% for some cases. 

Our calculation is performed on single PC (AMD 
Opteron(tm) Processor 248). Parallel implementation is 
easy for the current method. Since the major numeric 
cost is the computation of the matrix elements of the 
Hamiltonian during the search of the basis states, par- 
allel implementation can be simply realized by requiring 
each node handling some matrix elements of the Hamil- 
tonian. 

There are many possible ways to improve the current 
method. For example, it is worth to explore other type 
of basis states. In the present form, our method is an ex- 
tension to the mean field HF approximation. Mathemat- 
ically, the redundancy of the over-complete basis states is 
crucial for the convergence speed of the MP algorithm 0] . 
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By increase the redundancy of the over-complete basis 
states, i.e., enlarging the searching space, it is possible 
to speed up the convergence of MP method for searching 
the basis states. On the other hand, for particles moving 
in 3D space, storage of a single particle state needs siz- 
able memory, one may choose the basis states for single 
particles as product of one dimensional wave functions. 
In principle, this method is able to compute the exited 
states. With some modifications, it may be feasible to 
calculate the low-lying exited states. 



In summary, the current method is stable and free from 
the sign problem. It can apply to any system that can 
apply Hartree-Fock algorithm, and can be regarded as 
an extension to the Hartree-Fock algorithm. Several tens 
of determinant states are usually enough for meaningful 
result. This method may offer an alternative to explore 
quantum effects of Many body systems. 

This work is supported by the National Science Foun- 
dation of China (Grant No. 10375042). We acknowledge 
beneficial discussions with Prof. W. Wang. 
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